*****************************************
*** 2. LATE PSYCHOLOGICAL with CONTROLS
*** Prepare adjusted p-values
use "$gsdRawTemp\M1qvalue_ate_con_psy_1", clear
foreach num of numlist 2/5{
	append using "$gsdRawTemp\M1qvalue_ate_con_psy_`num'"
}
*append using "$gsdRawTemp\qvalue_atex_con_1"
*append using "$gsdRawTemp\qvalue_atex_con_2"
***Bonferroni q-values
qqvalue p , method(bonferroni) qvalue(qvalue)
gen bonferroni_q=round(qvalue, .001)
****Benjamini Hochbergen q-values
* Collect the total number of p-values tested
quietly sum p
local totalpvals = r(N)
* Sort the p-values in ascending order and generate a variable that codes each p-value's rank
quietly gen int original_sorting_order = _n
quietly sort p
quietly gen int rank = _n if p!=.
* Set the initial counter to 1 
local qval = 1
* Generate the variable that will contain the BH (1995) q-values
gen bh95_qval = 1 if p!=.
* Set up a loop that begins by checking which hypotheses are rejected at q = 1.000, then checks which hypotheses are rejected at q = 0.999, then checks which hypotheses are rejected at q = 0.998, etc.  The loop ends by checking which hypotheses are rejected at q = 0.001.
while `qval' > 0 {
	* Generate value qr/M
	quietly gen fdr_temp = `qval'*rank/`totalpvals'
	* Generate binary variable checking condition p(r) <= qr/M
	quietly gen reject_temp = (fdr_temp>=p) if fdr_temp~=.
	* Generate variable containing p-value ranks for all p-values that meet above condition
	quietly gen reject_rank = reject_temp*rank
	* Record the rank of the largest p-value that meets above condition
	quietly egen total_rejected = max(reject_rank)
	* A p-value has been rejected at level q if its rank is less than or equal to the rank of the max p-value that meets the above condition
	replace bh95_qval = `qval' if rank <= total_rejected & rank~=.
	* Reduce q by 0.001 and repeat loop
	quietly drop fdr_temp reject_temp reject_rank total_rejected
	local qval = `qval' - .001
}	
quietly sort original_sorting_order
*Save file and export to excel

preserve
keep if family == "behavioral_controls_t1"
mkmat bh95_qval, matrix(C3) rownames(parm)
restore
keep if family == "behavioral_controls_t2"
mkmat bh95_qval, matrix(C4) rownames(parm)


************************************
*** 3. LATE PSYCHOLOGICAL with GEOGRAPHY CONTROLS
*** Prepare adjusted p-values
use "$gsdRawTemp\M1qvalue_ate_con2_psy_1", clear
foreach num of numlist 2/5{
	append using "$gsdRawTemp\M1qvalue_ate_con2_psy_`num'"
}
*append using "$gsdRawTemp\qvalue_atex_geocon_1"
*append using "$gsdRawTemp\qvalue_atex_geocon_2"
***Bonferroni q-values
qqvalue p , method(bonferroni) qvalue(qvalue)
gen bonferroni_q=round(qvalue, .001)
****Benjamini Hochbergen q-values
* Collect the total number of p-values tested
quietly sum p
local totalpvals = r(N)
* Sort the p-values in ascending order and generate a variable that codes each p-value's rank
quietly gen int original_sorting_order = _n
quietly sort p
quietly gen int rank = _n if p!=.
* Set the initial counter to 1 
local qval = 1
* Generate the variable that will contain the BH (1995) q-values
gen bh95_qval = 1 if p!=.
* Set up a loop that begins by checking which hypotheses are rejected at q = 1.000, then checks which hypotheses are rejected at q = 0.999, then checks which hypotheses are rejected at q = 0.998, etc.  The loop ends by checking which hypotheses are rejected at q = 0.001.
while `qval' > 0 {
	* Generate value qr/M
	quietly gen fdr_temp = `qval'*rank/`totalpvals'
	* Generate binary variable checking condition p(r) <= qr/M
	quietly gen reject_temp = (fdr_temp>=p) if fdr_temp~=.
	* Generate variable containing p-value ranks for all p-values that meet above condition
	quietly gen reject_rank = reject_temp*rank
	* Record the rank of the largest p-value that meets above condition
	quietly egen total_rejected = max(reject_rank)
	* A p-value has been rejected at level q if its rank is less than or equal to the rank of the max p-value that meets the above condition
	replace bh95_qval = `qval' if rank <= total_rejected & rank~=.
	* Reduce q by 0.001 and repeat loop
	quietly drop fdr_temp reject_temp reject_rank total_rejected
	local qval = `qval' - .001
}	
quietly sort original_sorting_order
*Save file and export to excel

preserve
keep if family == "behavioral_controls2_t1"
mkmat bh95_qval, matrix(C5) rownames(parm)
restore
keep if family == "behavioral_controls2_t2"
mkmat bh95_qval, matrix(C6) rownames(parm)

*****************************************
*** 7. LATE SOCIOECONOMIC SURVEY Difference between T1 and T2
*** Prepare adjusted p-values
use "$gsdRawTemp\M1qvalue_ate_difference2_psy_1", clear
foreach num of numlist 2/5{
	append using "$gsdRawTemp\M1qvalue_ate_difference2_psy_`num'"
}
***Bonferroni q-values
qqvalue p , method(bonferroni) qvalue(qvalue)
gen bonferroni_q=round(qvalue, .001)
****Benjamini Hochbergen q-values
* Collect the total number of p-values tested
quietly sum p
local totalpvals = r(N)
* Sort the p-values in ascending order and generate a variable that codes each p-value's rank
quietly gen int original_sorting_order = _n
quietly sort p
quietly gen int rank = _n if p!=.
* Set the initial counter to 1 
local qval = 1
* Generate the variable that will contain the BH (1995) q-values
gen bh95_qval = 1 if p!=.
* Set up a loop that begins by checking which hypotheses are rejected at q = 1.000, then checks which hypotheses are rejected at q = 0.999, then checks which hypotheses are rejected at q = 0.998, etc.  The loop ends by checking which hypotheses are rejected at q = 0.001.
while `qval' > 0 {
	* Generate value qr/M
	quietly gen fdr_temp = `qval'*rank/`totalpvals'
	* Generate binary variable checking condition p(r) <= qr/M
	quietly gen reject_temp = (fdr_temp>=p) if fdr_temp~=.
	* Generate variable containing p-value ranks for all p-values that meet above condition
	quietly gen reject_rank = reject_temp*rank
	* Record the rank of the largest p-value that meets above condition
	quietly egen total_rejected = max(reject_rank)
	* A p-value has been rejected at level q if its rank is less than or equal to the rank of the max p-value that meets the above condition
	replace bh95_qval = `qval' if rank <= total_rejected & rank~=.
	* Reduce q by 0.001 and repeat loop
	quietly drop fdr_temp reject_temp reject_rank total_rejected
	local qval = `qval' - .001
}	
quietly sort original_sorting_order
*Save file and export to excel

mkmat bh95_qval, matrix(C8) rownames(parm)


*****************************************
*** 7. LATE SOCIOECONOMIC SURVEY Difference between T1 and T2
*** Prepare adjusted p-values
use "$gsdRawTemp\M1qvalue_ate_difference3_psy_1", clear
foreach num of numlist 2/5{
	append using "$gsdRawTemp\M1qvalue_ate_difference3_psy_`num'"
}
***Bonferroni q-values
qqvalue p , method(bonferroni) qvalue(qvalue)
gen bonferroni_q=round(qvalue, .001)
****Benjamini Hochbergen q-values
* Collect the total number of p-values tested
quietly sum p
local totalpvals = r(N)
* Sort the p-values in ascending order and generate a variable that codes each p-value's rank
quietly gen int original_sorting_order = _n
quietly sort p
quietly gen int rank = _n if p!=.
* Set the initial counter to 1 
local qval = 1
* Generate the variable that will contain the BH (1995) q-values
gen bh95_qval = 1 if p!=.
* Set up a loop that begins by checking which hypotheses are rejected at q = 1.000, then checks which hypotheses are rejected at q = 0.999, then checks which hypotheses are rejected at q = 0.998, etc.  The loop ends by checking which hypotheses are rejected at q = 0.001.
while `qval' > 0 {
	* Generate value qr/M
	quietly gen fdr_temp = `qval'*rank/`totalpvals'
	* Generate binary variable checking condition p(r) <= qr/M
	quietly gen reject_temp = (fdr_temp>=p) if fdr_temp~=.
	* Generate variable containing p-value ranks for all p-values that meet above condition
	quietly gen reject_rank = reject_temp*rank
	* Record the rank of the largest p-value that meets above condition
	quietly egen total_rejected = max(reject_rank)
	* A p-value has been rejected at level q if its rank is less than or equal to the rank of the max p-value that meets the above condition
	replace bh95_qval = `qval' if rank <= total_rejected & rank~=.
	* Reduce q by 0.001 and repeat loop
	quietly drop fdr_temp reject_temp reject_rank total_rejected
	local qval = `qval' - .001
}	
quietly sort original_sorting_order
*Save file and export to excel

mkmat bh95_qval, matrix(C9) rownames(parm)

